
import numpy as np

my_array = np.array([1, 2, 3, 4, 5]) # basic 1-D numpy array
print('my_array sum', my_array.sum())
print('my_array reshaped', my_array.reshape(1, 5)) # reshape into 2-D array

# numpy arrays can only contain one type of data
print('my_array dtype', my_array.dtype)

my_array sum 15
my_array reshaped [[1 2 3 4 5]]
my_array dtype int64


import pandas as pd

# basic pandas series - similar to a 1-D numpy array. series must contain data 
# of the same type. a series can be thought of as a single column of a table
my_series = pd.Series([10, 20, 30, 40, 50], name='my_series')
print('my_series sum', my_series.sum())
print('my_series mean', my_series.mean())

# create a pandas dataframe with some dummy data. DataFrames are comprised of multiple Series objects
my_df = pd.DataFrame({
    'county': ['Oktibbeha', 'Hinds', 'Rankin', 'Madison', 'Lauderdale'],
    'avg_income': [60000, 50000, 70000, 80000, 40000],
    'avg_age': [30, 35, 40, 45, 50],
    'population': [100000, 200000, 300000, 400000, 500000]
})

print('Oktibbeha county average age:', my_df[my_df['county'] == 'Oktibbeha']['avg_age'].values[0]) # select row containing data for Oktibbeha county and get the avg_age value
print(f'Total population across all {len(my_df)} counties:', my_df['population'].sum()) # sum the population column
print('Pearson correlation between avg_income and avg_age:', my_df['avg_income'].corr(my_df['avg_age'])) # calculate the correlation between avg_income and avg_age

my_series sum 150
my_series mean 30.0
Oktibbeha county average age: 30
Total population across all 5 counties: 1500000
Pearson correlation between avg_income and avg_age: -0.1


import matplotlib.pyplot as plt

# create a subplot with 1 row and 2 columns
fig, ax = plt.subplots(1, 2, figsize=(10, 5))

# plot the avg_income and avg_age columns on the first subplot
ax[0].scatter(my_df['avg_income'], my_df['avg_age'])
ax[0].set_xlabel('Average Income')
ax[0].set_ylabel('Average Age')
ax[0].set_title('Income vs Age')

# plot the avg_income and population columns on the second subplot
ax[1].scatter(my_df['avg_income'], my_df['population'])
ax[1].set_xlabel('Average Income')
ax[1].set_ylabel('Population')
ax[1].set_title('Income vs Population')

fig.tight_layout() # automatically adjust positions of figure elements to best fit everything
plt.show() # display the plot


import shapely.geometry

my_point = shapely.geometry.Point(-88.792972, 33.453323)
my_line = shapely.geometry.LineString([
    (-88.792972, 33.453333), 
    (-88.792982, 33.453343), 
    (-88.792992, 33.453333)
])

# polygon coordinates are same as line coordinates, but they form a closed shape (last point is same as first point)
my_polygon = shapely.geometry.Polygon([
    (-88.792972, 33.453333), 
    (-88.792982, 33.453343), 
    (-88.792992, 33.453333), 
    (-88.792972, 33.453333)
])

# shapely can help us calculate areas, lengths, and other geometric properties
print('Area of my_polygon:', my_polygon.area)
print('Length of my_line:', my_line.length)
print('Bounding box of my_polygon:', my_polygon.bounds)

# can also examine relationships between shapely objects
print('Does my_point intersect my_line?', my_point.intersects(my_line))
print('Distance between my_point and my_line:', my_point.distance(my_line))

# use ipython display funciton to show shapely objects in Jupyter notebooks
display(my_point)
display(my_line)
display(my_polygon)

# matplotlib can also be used to plot shapely objects
fig, ax = plt.subplots(3, 1, figsize=(10, 8))

ax[0].scatter(*my_point.xy, color='green')
ax[0].set_title('Shapely Point')

ax[1].plot(*my_line.xy, color='blue')
ax[1].set_title('Shapely Line')

ax[2].plot(*my_polygon.exterior.xy, color='red')
ax[2].set_title('Shapely Polygon')

for a in ax:
    a.set_aspect('equal') # make sure x and y axes have same scale
    a.set_xlabel('Longitude')
    a.set_ylabel('Latitude')

fig.tight_layout()
plt.show()

# polygons can also have interior holes
my_polygon_with_hole = shapely.geometry.Polygon(
    shell=[
        (-88.792972, 33.453333), 
        (-88.792982, 33.453343), 
        (-88.792992, 33.453333), 
        (-88.792972, 33.453333)
    ], 
    holes=[
        [
            (-88.792985, 33.453335), 
            (-88.792988, 33.453336),
            (-88.792987, 33.453337),
            (-88.792985, 33.453335),
        ], 
        [
            (-88.792982, 33.453335), 
            (-88.792984, 33.453336),
            (-88.792983, 33.453337),
            (-88.792982, 33.453335),
        ], 
    ]
)

fig, ax = plt.subplots(1, 1, figsize=(5, 5))

# working with polygons with holes can be a bit tricky, but matplotlib can plot them
ax.plot(*my_polygon_with_hole.exterior.xy, color='red') # plot the exterior first
for hole in my_polygon_with_hole.interiors:
    ax.plot(*hole.xy, color='blue') # not plot the holes

ax.set_aspect('equal')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_title('Shapely Polygon with Holes')

fig.tight_layout()
plt.show()

Area of my_polygon: 9.999999992137419e-11
Length of my_line: 2.828427123634253e-05
Bounding box of my_polygon: (-88.792992, 33.453333, -88.792972, 33.453343)
Does my_point intersect my_line? False
Distance between my_point and my_line: 1.0000000003174137e-05


from pyproj import Transformer
# you can optionally use the `pyproj.transform` function to convert between CRSs 
# without creating a `Transformer` object, though I typically use the latter 
# since the `Transformer` object can be reused for multiple conversions

# pyproj does not directly support shapely objects, but in recent version of 
# shapely a utility function was added to help with this
from shapely.ops import transform

transformer = Transformer.from_crs(
    'epsg:4326', # WGS 84 (i.e., standard lat/lon in degrees)
    'epsg:32616', # UTM Zone 16N (meters)
    always_xy=True, # longitude first, latitude second (already defined our geometries this way earlier)
    area_of_interest=None, # can optioanally define an area of interest, not covered here
)

my_point_transformed = transform(transformer.transform, my_point)
my_line_transformed = transform(transformer.transform, my_line)
my_polygon_transformed = transform(transformer.transform, my_polygon)

# now we can look at properties of the transformed objects in meters
print('Area of my_polygon_transformed:', my_polygon_transformed.area, 'meters')
print('Length of my_line_transformed:', my_line_transformed.length, 'meters')
print('Distance between my_point_transformed and my_line_transformed:', my_point_transformed.distance(my_line_transformed), 'meters')

# now we can plot the transformed objects - note that the x and y axes are different
fig, ax = plt.subplots(3, 1, figsize=(10, 8))

ax[0].scatter(*my_point_transformed.xy, color='green')
ax[0].set_title('Shapely Transformed Point')

ax[1].plot(*my_line_transformed.xy, color='blue')
ax[1].set_title('Shapely Transformed Line')

ax[2].plot(*my_polygon_transformed.exterior.xy, color='red')
ax[2].set_title('Shapely Transformed Polygon')

for a in ax:
    a.set_aspect('equal') # make sure x and y axes have same scale
    a.set_xlabel('X (meters)')
    a.set_ylabel('Y (meters)')

fig.tight_layout()
plt.show()

Area of my_polygon_transformed: 1.031061426538567 meters
Length of my_line_transformed: 2.8943421527710393 meters
Distance between my_point_transformed and my_line_transformed: 1.109061591448836 meters


import geopandas as gpd

my_gs = gpd.GeoSeries([my_point, my_line, my_polygon], crs='4326') # make sure we define the CRS

# union of all geometries in the GeoSeries
my_gs_union = my_gs.unary_union
display(my_gs_union)

# Like Pandas, GeoPandas plotting functions that can be used to visualize data
gpd.plotting.plot_series(my_gs, edgecolor='black', facecolor='none', linewidth=2)

# we can create a GeoDataFrame to assign attributes to each geometry
my_gdf = gpd.GeoDataFrame(
    data={
        'name': ['my_point', 'my_line', 'my_polygon'],
        'value': [10, 20, 30]
    },
    geometry=my_gs
)

# Now we can create a plot of the GeoDataFrame and color the geometries based on the 'value' column
my_gdf.plot(column='value', legend=True)
plt.show()

# geopandas can load data from shapefiles
counties_gdf = gpd.read_file('data/cb_2018_us_county_20m/cb_2018_us_county_20m.shp')

display(counties_gdf.head())
# we can manipulate the data in the GeoDataFrame like we would a Pandas DataFrame

ms_counties_gdf = counties_gdf[counties_gdf['STATEFP'] == '28'] # select only Mississippi counties (FP code 28)
ms_counties_gdf.plot()

# Geopandas can convert between CRSs using the `to_crs` method
ms_counties_gdf = ms_counties_gdf.to_crs('epsg:32616') # convert to UTM Zone 16N

# plot to show the difference in the x and y axes
ms_counties_gdf.plot()
plt.show()

# We can join GeoDataFrames to other DataFrames using the `merge` method
my_df['county'] = my_df['county'].str.title() # make sure the county names are capitalized
ms_county_data = ms_counties_gdf.merge(
    my_df, # DataFrame to merge
    left_on='NAME', # column in left GeoDataFrame to join on
    right_on='county', # column in right DataFrame to join on
    how='inner' # type of join to perform
)
display(ms_county_data.head())

# Now we can make a fancy plot of the data
fig, ax = plt.subplots(1, 1, figsize=(10, 10))

# plot the Mississippi counties, pass in the `ax` object so we can overlay other 
# plots and configure the plot
ms_counties_gdf.plot(ax=ax, facecolor='none',) # plot all counties
ms_county_data.plot(ax=ax, column='population') # plot counties we have data for, color by population

ax.set_aspect('equal') # make sure x and y axes have same scale
ax.set_xlabel('X (meters)')
ax.set_ylabel('Y (meters)')
ax.set_title('Population of Mississippi Counties')

fig.tight_layout()
plt.show()

# Can look at spatial statistics as well
print('Average area of Mississippi counties:', ms_counties_gdf.area.mean(), 'meters squared')
print('Total area of Mississippi counties:', ms_counties_gdf.area.sum(), 'meters squared')

Average area of Mississippi counties: 1505025316.27986 meters squared
Total area of Mississippi counties: 123412075934.94852 meters squared


from osgeo import gdal

# open a raster dataset
dataset = gdal.Open('./data/m_3308834_ne_16_060_20230820.jp2')

# look at dataset properties

print(f'Size is {dataset.RasterXSize} x {dataset.RasterYSize} x {dataset.RasterCount}')
print(f'Raster has {dataset.RasterCount} bands') # Number of bands in the raster
print(f'Projection is {dataset.GetProjection()}') # CRS of the raster
print(f'Geotransform is {dataset.GetGeoTransform()}') # Will go over this later

# read the data from the raster into a numpy array
band = dataset.GetRasterBand(1).ReadAsArray() # get the first band (red in this case

print(f'Shape of the band array is {band.shape}')
print(f'Band data type is {band.dtype}')
print(f'Band min-max is {band.min()} - {band.max()}')

fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.imshow(band, cmap='gray') # plot the data
plt.show()

# read all bands into a 3-D numpy array
data = np.stack([dataset.GetRasterBand(i).ReadAsArray() for i in range(1, dataset.RasterCount + 1)], axis=-1) # iterate over all bands and stack them along the last axis

print(f'Shape of the data array is {data.shape}')
print(f'Data type is {data.dtype}')
print(f'Data min-max is {data.min()} - {data.max()}')

# plot the data
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.imshow(data) # plot the data
plt.show()

Size is 11599 x 13904 x 4
Raster has 4 bands
Projection is PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]
Geotransform is (-9886563.026870845, 0.6000000000000321, 0.0, 3961861.2086821054, 0.0, -0.5999999999999933)
Shape of the band array is (13904, 11599)
Band data type is uint8
Band min-max is 0 - 255

Shape of the data array is (13904, 11599, 4)
Data type is uint8
Data min-max is 0 - 255


# load the affine transformtation matrix

transform = dataset.GetGeoTransform()
print(f'Geotransform is {transform}')
print(f'Affine transformation matrix in 3x3 form is\n{np.array([[transform[1], transform[2], transform[0]], [transform[4], transform[5], transform[3]], [0, 0, 1]])}')
print(f'Pixel size is {transform[1]} x {transform[5]} meters') # pixel size in y direction is negative by convention

# close the dataset - we are done for now
dataset = None

Geotransform is (-9886563.026870845, 0.6000000000000321, 0.0, 3961861.2086821054, 0.0, -0.5999999999999933)
Affine transformation matrix in 3x3 form is
[[ 6.00000000e-01  0.00000000e+00 -9.88656303e+06]
 [ 0.00000000e+00 -6.00000000e-01  3.96186121e+06]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
Pixel size is 0.6000000000000321 x -0.5999999999999933 meters


import rasterio as rio
from rasterio.plot import show # rasterio convenience function for plotting rasters

# use context manager to open a raster dataset - similar to standard Python file I/O
with rio.open('./data/m_3308834_ne_16_060_20230820.jp2') as dataset:
    data = dataset.read() # read all bands into a 3-D numpy array
    nir_band = dataset.read(4) # read the 4th band (NIR in this case)
    meta = dataset.meta # get metadata (CRS, transform, dtype, shape, etc.)
    transform = dataset.transform # get the affine transformation matrix directly (can also use `meta['transform']`)

print(f'Shape of the data array is {data.shape}')
print(f'Data type is {data.dtype}')
print(f'Data min-max is {data.min()} - {data.max()}')
print('Data metadata:', meta)
print('Affine transformation matrix:', transform)

# rio contains a convenience function for plotting rasters
fig, ax = plt.subplots(1, 1, figsize=(10, 10))

show(nir_band, ax=ax, cmap='viridis', transform=transform) # plot the NIR band

fig.suptitle('NAIP - NIR reflectance for MSU and Oktibehhah county', fontsize=16)
plt.show()

# show a false color composite of the data
fig, ax = plt.subplots(1, 1, figsize=(10, 10))

show(data[[3, 0, 1], ...], ax=ax, transform=transform) # plot the NIR, R, and G bands
fig.suptitle('NAIP - False color composite for MSU and Oktibehhah county', fontsize=16)
plt.show()

Shape of the data array is (4, 13904, 11599)
Data type is uint8
Data min-max is 0 - 255
Data metadata: {'driver': 'JP2OpenJPEG', 'dtype': 'uint8', 'nodata': None, 'width': 11599, 'height': 13904, 'count': 4, 'crs': CRS.from_epsg(3857), 'transform': Affine(0.6000000000000321, 0.0, -9886563.026870845,
       0.0, -0.5999999999999933, 3961861.2086821054)}
Affine transformation matrix: | 0.60, 0.00,-9886563.03|
| 0.00,-0.60, 3961861.21|
| 0.00, 0.00, 1.00|


water_bodies_gdf = gpd.read_file('./data/water_bodies.gpkg')

# IMPORTANT - remember we need to make sure the CRS of the GeoDataFrame matches the CRS of the raster
water_bodies_gdf = water_bodies_gdf.to_crs(meta['crs'])

fig, ax = plt.subplots(1, 1, figsize=(10, 10))

# plot the raster first in color infrared
show(data[[3, 0, 1], ...], ax=ax, transform=transform) # plot the NIR, R, and G bands
water_bodies_gdf.plot(ax=ax, facecolor='blue', edgecolor='black') # overlay the water bodies

fig.suptitle('NAIP - False color composite with water bodies', fontsize=16)
fig.tight_layout()
plt.show()


# partition the raster into 256x256 pixel tiles - we can do this in Numpy
import os

tile_size = 256
data_output_dir = './data/tiles/grid/input'
os.makedirs(data_output_dir, exist_ok=True)

with rio.open('./data/m_3308834_ne_16_060_20230820.jp2') as dataset:
    for x_offset in range(0, dataset.width, tile_size):
        for y_offset in range(0, dataset.height, tile_size):
            
            window = rio.windows.Window(
                col_off=x_offset,
                row_off=y_offset,
                width=tile_size,
                height=tile_size,
            )
            
            # only read in NIR, Red, and Green bands
            tile = dataset.read((4, 1, 2), window=window)
            
            if tile.shape[1:] != (tile_size, tile_size):
                continue # skip incomplete tiles
            
            if dataset.nodata in tile:
                continue # skip tiles with nodata values
            
            if tile.min() == tile.max():
                continue # skip tiles with no variation
            
            # update the metadata to reflect the properties of the new tile
            new_meta = meta.copy()
            new_meta['transform'] = rio.windows.transform(
                window,  # provide the window to get the transform
                dataset.transform # transform of the original raster
            ) # returns the transform of the new window
            new_meta['count'], new_meta['height'], new_meta['width'] = tile.shape # update the height and width
            new_meta['driver'] = 'GTiff' # update the driver to GeoTIFF
            new_meta['compress'] = 'lzw' # use LZW compression to minimize file size

            # save the raster 
            with rio.open(f'{data_output_dir}/tile_{x_offset}_{y_offset}.tif', 'w', **new_meta) as dst:
                dst.write(tile)


from rasterio.features import rasterize
from glob import glob

target_dir = './data/tiles/grid/output'
os.makedirs(target_dir, exist_ok=True)

input_tiles = glob(f'{data_output_dir}/*.tif')

# keep track of tiles with data
tiles_with_data = []
for tile in input_tiles:
    
    # need to extract the x and y offsets from the tile filename
    tile_file_base = os.path.basename(tile).split('.')[0] # use .split to remove the file extension
    x_offset = tile_file_base.split('_')[1]
    y_offset = tile_file_base.split('_')[2]
    
    with rio.open(tile) as src:
        tile_bounds = src.bounds # get the bounds of the tile
    
    # print(*tile_bounds)
    
    tile_features = water_bodies_gdf.clip(tile_bounds) # clip the water bodies to the tile bounds
    
    out_file = f'{target_dir}/tile_{x_offset}_{y_offset}_water_bodies.tif'
    
    # handle the case where there are no water bodies in the tile
    if tile_features.empty:
        # create an empty numpy array with the same shape as the tile
        target_arr = np.zeros((src.height, src.width), dtype=np.uint8) # use uint8 to save space (don't need 32-bit integers)
        # print(target_arr.shape)
    
    else: # rasterize the water bodies
        target_arr = rasterize(
            [(geom, 1) for geom in tile_features.geometry], # list of geometries to rasterize
            out_shape=(src.height, src.width), # shape of the output array
            transform=src.transform, # affine transformation matrix
            fill=0, # fill value for pixels not covered by geometries
            all_touched=True, # use all pixels touched by the geometry
            default_value=1, # value to assign to pixels covered by geometries
            dtype=np.uint8 # data type of the output array
        )
        # print(target_arr.shape)
        # print(target_arr)
        tiles_with_data.append(out_file)
    
    # update the metadata to for the new target raster
    new_meta = src.meta.copy()
    new_meta['count'] = 1 # only one band
    new_meta['dtype'] = 'uint8' # data type is unsigned 8-bit integer
    new_meta['compress'] = 'lzw' # use LZW compression to minimize file size
    
    # save the raster
    with rio.open(out_file, 'w', **new_meta) as dst:
        dst.write(target_arr, 1) # write the data to the first band


# lets visualize the output
out_tile = tiles_with_data[0] # get the first tile with data
x_offset = out_tile.split('_')[1] # extract the x and y offsets from the filename to load the corresponding input tile
y_offset = out_tile.split('_')[2]
in_tile = f'./data/tiles/grid/input/tile_{x_offset}_{y_offset}.tif'

# load the input and output tiles
with rio.open(in_tile) as src:
    in_data = src.read()
    in_transform = src.transform
    
with rio.open(out_tile) as src:
    out_data = src.read()
    out_transform = src.transform

# create a data mask to overlay the water bodies
data_mask = np.ma.masked_where(out_data == 0, out_data) # mask the data where the water bodies are not present

fig, ax = plt.subplots(1, 1, figsize=(5, 5))

show(in_data, ax=ax,) # plot the input imagery
ax.imshow(data_mask[0], cmap='Blues_r') # overlay the water bodies

plt.show()


bbox_out_dir = './data/tiles/grid/bbox'
os.makedirs(bbox_out_dir, exist_ok=True)

for tile_file in tile_files:
    
    with rio.open(tile_file) as src:
        tile_bounds = src.bounds # get the bounds of the tile
        data = src.read() # read the data
        
        tile_size = data.shape # get the tile size
    
    # clip the water bodies to the tile bounds
    tile_features = water_bodies_gdf.clip(tile_bounds)
    # print(tile_features)
    
    if tile_features.empty:
        continue
    
    # assign a class label to the tile features - 0 for now since we only have water bodies
    tile_features['class'] = 0
    
    # for each water body in the tile, extract the x_center and y_center and scale it [0, 1] with respect to the tile
    for i, row in tile_features.iterrows():
        
        left_bound = row.geometry.bounds[0] # minx
        right_bound = row.geometry.bounds[2] # maxx
        bottom_bound = row.geometry.bounds[1] # miny
        top_bound = row.geometry.bounds[3] # maxy
        
        # scale bounds to [0, 1] with respect to the tile
        tile_width = tile_bounds.right - tile_bounds.left
        tile_height = tile_bounds.top - tile_bounds.bottom
        
        left_bound = (left_bound - tile_bounds.left) / tile_width
        right_bound = (right_bound - tile_bounds.left) / tile_width
        top_bound = (top_bound - tile_bounds.bottom) / tile_height
        bottom_bound = (bottom_bound - tile_bounds.bottom) / tile_height
        
        # now determine center of the bounding box
        x_center = (left_bound + right_bound) / 2
        y_center = (bottom_bound + top_bound) / 2
        
        # get the width and height of the bounding box
        width = right_bound - left_bound
        height = top_bound - bottom_bound
        
        # lets visualize the resultant bounding box
        fig, ax = plt.subplots(1, 1, figsize=(5, 5))
        show(data, ax=ax) # plot the input imagery
        # overlay bounding box
        ax.add_patch(
            plt.Rectangle(
                (
                    (x_center - (width / 2)) * tile_size[1] - 1,
                    tile_size[2] - ((y_center + (height / 2)) * tile_size[2]), # flip the y-axis (to match the rasterio plot convention)
                ), 
                width * tile_size[1], 
                height * tile_size[2], 
                edgecolor='red', 
                facecolor='none',
            )
        )
        
        # overlay annotations just to make sure we are getting the right bounding box
        with rio.open(tile_file.replace('input', 'output').replace('.tif', '_water_bodies.tif')) as src:
            show(src.read(), ax=ax, cmap='Blues', alpha=0.5)
        
        plt.show() # NOTE: bounding box is not perfect in visualization - but is correct in calculations
        
        # now we can save the bounding box information to a txt file in the format
        # label_idx x_center y_center width height
        print(f'{row["class"]} {x_center} {y_center} {width} {height}')
        with open(f'{bbox_out_dir}/{os.path.basename(tile_file).replace(".tif", ".txt")}', 'w') as f:
            f.write(f'{row["class"]} {x_center} {y_center} {width} {height}')

0 0.7125272063136132 0.14650894787631447 0.4982814288462839 0.29301789575262893

0 0.602858935958412 0.2897105415981728 0.7942821280831759 0.5794210831963456

0 0.3286553619675121 0.3823663766462036 0.41581593102313685 0.18237540834877874

0 0.03837881125955604 0.7406687182442567 0.07675762251911208 0.08847924761837422

0 0.1484215167090307 0.04313378509775361 0.26060225033448803 0.08626757019550722

0 0.03029451550302303 0.008078537460576579 0.06058903100604606 0.016157074921153158

0 0.5 0.5 1.0 1.0

0 0.02335210438110359 0.6883277670942405 0.04670420876220718 0.4212775367732797

0 0.5 0.5 1.0 1.0

0 0.9048944775674921 0.21040464956652286 0.19021104486501572 0.4208092991330457

	STATEFP	COUNTYFP	COUNTYNS	AFFGEOID	GEOID	NAME	LSAD	ALAND	AWATER	geometry
0	37	017	01026336	0500000US37017	37017	Bladen	06	2265887723	33010866	POLYGON ((-78.90200 34.83527, -78.79960 34.850...
1	37	167	01025844	0500000US37167	37167	Stanly	06	1023370459	25242751	POLYGON ((-80.49737 35.20210, -80.29542 35.502...
2	39	153	01074088	0500000US39153	39153	Summit	06	1069181981	18958267	POLYGON ((-81.68699 41.13596, -81.68495 41.277...
3	42	113	01213687	0500000US42113	42113	Sullivan	06	1165338428	6617028	POLYGON ((-76.81373 41.59003, -76.22014 41.541...
4	48	459	01384015	0500000US48459	48459	Upshur	06	1509910100	24878888	POLYGON ((-95.15274 32.66095, -95.15211 32.902...

Working with Geospatial Data in Python¶

Background¶

Common libraries for data and numerical analysis¶

Numpy¶

Pandas¶

Matplotlib¶

Geospatial libraries¶

Vector data¶

Shapely¶

A note on geometry coordinates¶

Pyproj¶

Geographic Coordinate Systems and their limitations¶

A note on GCSs¶

Projected Coordinate Systems¶

Coordinate Reference Systems¶

Converting between CRSs using Pyproj¶

Geopandas¶

Raster data¶

GDAL¶

Affine transforms¶

Rasterio¶

Example - Clipping a raster and vector data for deep learning training¶

First - loading the vector data¶

Clipping the data by grids¶

Extracting bounding boxes for object detection¶

Conclusion¶

	STATEFP	COUNTYFP	COUNTYNS	AFFGEOID	GEOID	NAME	LSAD	ALAND	AWATER	geometry	county	avg_income	avg_age	population
0	28	049	00695748	0500000US28049	28049	Hinds	06	2252697076	19450105	POLYGON ((148567.197 3571491.877, 151245.417 3...	Hinds	50000	35	200000
1	28	105	00695775	0500000US28105	28105	Oktibbeha	06	1186831624	9568768	POLYGON ((306070.347 3712277.418, 306071.132 3...	Oktibbeha	60000	30	100000
2	28	089	00695768	0500000US28089	28089	Madison	06	1850581711	71295137	POLYGON ((176088.403 3609292.360, 180333.943 3...	Madison	80000	45	400000
3	28	121	00695781	0500000US28121	28121	Rankin	06	2008615552	79136591	POLYGON ((193000.253 3561837.331, 198777.758 3...	Rankin	70000	40	300000
4	28	075	00695761	0500000US28075	28075	Lauderdale	06	1822458570	29780252	POLYGON ((320291.112 3606006.706, 369602.675 3...	Lauderdale	40000	50	500000